Tutorial Part VI — Population analysis in R
Introduction
This is a brief tutorial about constructing and analyzing relatively simple Cormack-Jolly-Seber capture-mark-recapture (CMR) models in R. For this we use the R package “marked” (Laake et al. 2013).
First we get an overview of the data, then we develop different models and finally we analyse and visualize the results.
Setup
Load (and install) the following packages:
package.list=c("here",
"marked",
"skimr",
"tidyverse",
"sf",
"tmap"
)
for (package in package.list) {
if (!require(package, character.only=T, quietly=T)) {
install.packages(package)
library(package, character.only=T)
}
}set the theme for some ggplots:
theme_set(
theme(
panel.grid = element_blank(),
panel.background = element_rect(fill = "white"),
panel.border = element_rect(
colour = "black",
fill = NA,
linewidth = 1
),
axis.ticks = element_line(colour = "black", linewidth = 1),
axis.ticks.length = unit(0.2, "cm")
)
)Load and manipulate data
Here, we want to directly set some columns as factors for all subsequent analyses:
starlings <-
read_csv(here("data", "data_berlin", "animal_data", "starlings.csv")) %>%
mutate(
sex = as.factor(sex),
habitat = as.factor(habitat),
infected = as.factor(infected)
) %>%
rename(ID = "...1")
# Short overview
skim(starlings)| Name | starlings |
| Number of rows | 605 |
| Number of columns | 5 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| factor | 3 |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| ch | 0 | 1 | 7 | 7 | 0 | 32 | 0 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| sex | 0 | 1.00 | FALSE | 2 | Fem: 323, Mal: 282 |
| habitat | 0 | 1.00 | FALSE | 2 | urb: 311, rur: 294 |
| infected | 202 | 0.67 | FALSE | 2 | Yes: 214, No: 189 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| ID | 0 | 1 | 303 | 174.79 | 1 | 152 | 303 | 454 | 605 | ▇▇▇▇▇ |
set.seed(123)Each row in this table represents one individual. We have 5 columns in this dataset: ID, the individual identifier, ch, the capture history (0 = undetected, 1 = detected),the sex of the captured individuals, the habitat type the birds were captured, and health status, i.e. if birds were infected
Explore capture locations
#load environmental data on capturing
locations <-
read_csv(here(
"data",
"data_berlin",
"animal_data",
"starlings_capture_location_data.csv"
))
# Short overview
skim(locations)| Name | locations |
| Number of rows | 2 |
| Number of columns | 10 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| numeric | 9 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| Habitat | 0 | 1 | 5 | 5 | 0 | 2 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Lon | 0 | 1 | 13.57 | 0.19 | 13.44 | 13.51 | 13.57 | 13.64 | 13.71 | ▇▁▁▁▇ |
| Lat | 0 | 1 | 52.44 | 0.05 | 52.40 | 52.42 | 52.44 | 52.45 | 52.47 | ▇▁▁▁▇ |
| temperature_2016 | 0 | 1 | 30.25 | 1.06 | 29.50 | 29.88 | 30.25 | 30.62 | 31.00 | ▇▁▁▁▇ |
| temperature_2017 | 0 | 1 | 39.00 | 5.66 | 35.00 | 37.00 | 39.00 | 41.00 | 43.00 | ▇▁▁▁▇ |
| temperature_2018 | 0 | 1 | 33.50 | 3.54 | 31.00 | 32.25 | 33.50 | 34.75 | 36.00 | ▇▁▁▁▇ |
| temperature_2019 | 0 | 1 | 28.50 | 2.12 | 27.00 | 27.75 | 28.50 | 29.25 | 30.00 | ▇▁▁▁▇ |
| temperature_2020 | 0 | 1 | 28.50 | 0.71 | 28.00 | 28.25 | 28.50 | 28.75 | 29.00 | ▇▁▁▁▇ |
| temperature_2021 | 0 | 1 | 30.00 | 1.41 | 29.00 | 29.50 | 30.00 | 30.50 | 31.00 | ▇▁▁▁▇ |
| temperature_2022 | 0 | 1 | 29.50 | 1.41 | 28.50 | 29.00 | 29.50 | 30.00 | 30.50 | ▇▁▁▁▇ |
loc_sf <- st_as_sf(locations, coords = c("Lon", "Lat"), crs = 4326)
tmap_mode("view")
tm_shape(loc_sf) + tm_dots(col = "Habitat", size = 2) Plot the capturing data
# Capture histories
starlings %>%
ggplot(mapping = aes(x = ch)) +
geom_bar(colour = "#262674", fill = "#262674") +
scale_y_continuous(
limits = c(0, 40),
breaks = seq(0, 40, by = 10),
expand = c(.001, .001)
) +
labs(x = "capture histories") +
theme(axis.text.x = element_text(angle = 45))# Sex
starlings %>%
ggplot(mapping = aes(x = sex, fill = sex)) +
geom_bar() +
scale_fill_manual(values = c(Female = "salmon", Male = "lightblue"))# Sex and capture histories
starlings %>%
ggplot(mapping = aes(x = ch, fill = sex)) +
geom_bar(position = "stack") +
scale_y_continuous(
limits = c(0, 40),
breaks = seq(0, 40, by = 10),
expand = c(.001, .001)
) +
labs(x = "capture histories", fill = "sex") +
theme(axis.text.x = element_text(angle = 45)) +
scale_fill_manual(values = c(Female = "salmon", Male = "lightblue"))Plot environmental data
locations_long <- locations %>%
pivot_longer(names_to = "year",
values_to = "temp",
cols = starts_with("Temp"))
ggplot(locations_long, aes(x = year, y = temp, group = Habitat)) +
geom_line(aes(col = Habitat), linewidth = 2) +
scale_colour_manual(values = c(urban = "darkgoldenrod", rural = "darkblue"))Basic Cormack-Jolly-Seber models
The default model uses constant detection and survival probabilities and constant time intervals between captures. We fit the model and have a look at the results. Please note: The estimates are on a logit scale
#models require dataframes, tibbles produce errors
starlings <- starlings %>%
as.data.frame()
cjs_m1 <- crm(data = starlings)
cjs_m1##
## crm Model Summary
##
## Npar : 2
## -2lnL: 1291.825
## AIC : 1295.825
##
## Beta
## Estimate
## Phi.(Intercept) 0.06757844
## p.(Intercept) 2.22899510
Npar = number of parameters -2lnL = AIC = Akaikes Information Criterion Phi = Survival probability between capture events p = Capturing probability per capture event
Estimate parameters on real scale
Write a function to assess results more intuitively on the real scale instead of the logit-scale. We can just paste our results from the previous model output
inverse_logit <- function(x) {
exp(x) / (1 + exp(x))
}
#e.g.
inverse_logit(0.06757844) # this is our survival probability## [1] 0.5168882
inverse_logit(2.22899510) # this is our survival probability## [1] 0.9028232
# alternative way to write this
inverse_logit(cjs_m1$results$beta$Phi)## (Intercept)
## 0.5168882
inverse_logit(cjs_m1$results$beta$p)## (Intercept)
## 0.9028232
Refit the model with confidence intervals
cjs_m1_2 <- cjs.hessian(cjs_m1)
cjs_m1_2##
## crm Model Summary
##
## Npar : 2
## -2lnL: 1291.825
## AIC : 1295.825
##
## Beta
## Estimate se lcl ucl
## Phi.(Intercept) 0.06757844 0.07249005 -0.07450207 0.2096589
## p.(Intercept) 2.22899510 0.24951951 1.73993685 2.7180533
Predict
Now we transform the logit estimates to real estimates by using the inverse logit
predict(object = cjs_m1_2,
newdata = data.frame(sex = c("Female", "Male")),
se = TRUE)## $Phi
## occ estimate se lcl ucl
## 1 1 0.5168882 0.01810184 0.4813831 0.5522236
##
## $p
## occ estimate se lcl ucl
## 1 2 0.9028232 0.02189121 0.850679 0.9380836
What does that mean? The survival probability between capture events was 0.52. The probability to capture an individual during a capture event was 0.90.
As mentioned before, we assume constant time intervals between capturing events. In reality, this is often not the case, e.g. capturing events have to be postponed due to weather conditions. Now we build a model where we take this into account.
CJS with irregular sampling intervals
cjs_m1_time <- crm(data = starlings,
# time.interval vector is the interval between each capture event.
# for 7 capture events, there are 6 intervals
time.intervals = c(1, 1, 1, 1.6, 2, 3))
predict(object = cjs_m1_time)## $Phi
## occ estimate
## 1 6 0.688108
##
## $p
## occ estimate
## 1 7 0.8858232
Now we have a higher survival probability of 0.69 between capturing events and a lower capturing probability per capturing event of 0.89
Fit models for model selection
Normally, there are different hypotheses about the factors that influence the survival probability of individuals. These can be characteristics of the individuals, such as sex, age or weight, or characteristics of the environment of the individuals, such as sealing, temperature or the number of nesting opportunities. For this purpose, different models are built that take the respective factors (explanatory or predictor variables) into account.
For a flexible approach we process the data, create design data to add variables that depend on time and fit the different models.
Group data
We process the data and set the grouping variables
m2_proc <- process.data(data = starlings,groups = c("sex","habitat","infected"))Create the design data
m2_design <- make.design.data(data = m2_proc)Model formulas
Build formulas for each parameter
phi_dot <- list(formula = ~ 1)
phi_sex <- list(formula = ~ sex)
p_sex <- list(formula = ~ sex)Fit a model based on the design data with constant survival and different detection probabilities between sexes.
cjs_m2 <- crm(
m2_proc,
m2_design,
model.parameters = list(Phi = phi_dot,
p = p_sex),
accumulate = FALSE
)##
Number of evaluations: 100 -2lnl: 1289.126943
cjs_m2##
## crm Model Summary
##
## Npar : 3
## -2lnL: 1289.117
## AIC : 1295.117
##
## Beta
## Estimate
## Phi.(Intercept) 0.0793894
## p.(Intercept) 1.7877066
## p.sexMale 0.7906268
predict(object = cjs_m2)## $Phi
## occ estimate
## 1 6 0.5198369
##
## $p
## sex occ estimate
## 1 Female 7 0.8566459
## 2 Male 7 0.9294541
Fit another model based on the design data with different survival and detection probabilities between sexes
cjs_m3 <- crm(
m2_proc,
m2_design,
model.parameters = list(Phi = phi_sex,
p = p_sex),
accumulate = FALSE
)##
Number of evaluations: 100 -2lnl: 1283.680512
Number of evaluations: 200 -2lnl: 1283.680512
cjs_m3##
## crm Model Summary
##
## Npar : 4
## -2lnL: 1283.681
## AIC : 1291.681
##
## Beta
## Estimate
## Phi.(Intercept) -0.09718598
## Phi.sexMale 0.34094476
## p.(Intercept) 1.98842947
## p.sexMale 0.46984162
predict(object = cjs_m3)## $Phi
## sex occ estimate
## 1 Female 6 0.4757226
## 2 Male 6 0.5606397
##
## $p
## sex occ estimate
## 1 Female 7 0.8795769
## 2 Male 7 0.9211642
Fit another model based on the design data with time-dependent survival and detection probability
cjs_m4 <- crm(
m2_proc,
m2_design,
model.parameters = list(Phi = list(formula = ~ time),
p = list(formula = ~ time)),
accumulate = FALSE
)##
Number of evaluations: 100 -2lnl: 1244.670707
Number of evaluations: 200 -2lnl: 1239.725494
Number of evaluations: 300 -2lnl: 1239.319519
Number of evaluations: 400 -2lnl: 1239.303701
Number of evaluations: 500 -2lnl: 1239.350357
Number of evaluations: 600 -2lnl: 1239.417367
Number of evaluations: 700 -2lnl: 1243.575689
Number of evaluations: 800 -2lnl: 1239.311438
Number of evaluations: 900 -2lnl: 1241.176613
Number of evaluations: 1000 -2lnl: 1239.329909
Number of evaluations: 1100 -2lnl: 1239.303675
cjs_m4##
## crm Model Summary
##
## Npar : 12
## -2lnL: 1239.304
## AIC : 1263.304
##
## Beta
## Estimate
## Phi.(Intercept) 1.5987971
## Phi.time2 -2.6225405
## Phi.time3 -1.8076489
## Phi.time4 -1.1463136
## Phi.time5 -1.3796043
## Phi.time6 -0.6965196
## p.(Intercept) 1.2003709
## p.time3 1.5729892
## p.time4 1.2288969
## p.time5 0.7978441
## p.time6 1.2524201
## p.time7 -0.2257742
predict(object = cjs_m4)## $Phi
## time occ estimate
## 1 6 6 0.7114173
## 2 5 5 0.5545798
## 3 4 4 0.6112295
## 4 3 3 0.4479760
## 5 2 2 0.2642989
## 6 1 1 0.8318502
##
## $p
## time occ estimate
## 1 7 7 0.7260348
## 2 6 6 0.9207653
## 3 5 5 0.8806095
## 4 4 4 0.9190321
## 5 3 3 0.9412192
## 6 2 2 0.7685908
Multiple models + model selection (AIC)
When testing multiple hypotheses, you can speed up model formulation for multiple models. Caution! Models should always be biologically meaningful, don’t just fit too many models without thinking about them.
fit.starlings.cjs.models <- function() {
# Apparent survival (Phi) formula
Phi.time <- list(formula = ~ time)
Phi.sex <- list(formula = ~ sex)
Phi.sex.time <-
list(formula = ~ sex * time) # interaction of sex and time
Phi.habitat.time <- list(formula = ~ habitat * time)
Phi.dot <- list(formula = ~ 1) # constant survival
# Detection probability (p) formula
p.sex <- list(formula = ~ sex) # differs between males and females
p.time <-
list(formula = ~ time) # one discrete estimate of p per capture event
p.habitat <- list(formula = ~ habitat)
p.dot <- list(formula = ~ 1) # constant detection
# Construct all combinations and put into one model table
cml <-
create.model.list(c("Phi", "p")) # makes all possibile combinations of those parameter formulas
results <- crm.wrapper(
cml,
data = m2_proc,
ddl = m2_design,
external = FALSE,
accumulate = FALSE,
hessian = TRUE
)
return(results)
}
# Run function
starlings_cjs_models <- fit.starlings.cjs.models()## Phi.dot.p.dot
## Phi.dot.p.habitat
##
Number of evaluations: 100 -2lnl: 1291.647716Phi.dot.p.sex
##
Number of evaluations: 100 -2lnl: 1289.126943Phi.dot.p.time
##
Number of evaluations: 100 -2lnl: 1286.513634
Number of evaluations: 200 -2lnl: 1284.979279
Number of evaluations: 300 -2lnl: 1284.814101
Number of evaluations: 400 -2lnl: 1284.787142
Number of evaluations: 500 -2lnl: 1284.8424
Number of evaluations: 600 -2lnl: 1285.724993
Number of evaluations: 700 -2lnl: 1284.7819
Number of evaluations: 100 -2lnl: 1284.88038
Number of evaluations: 200 -2lnl: 1284.787694Phi.habitat.time.p.dot
##
Number of evaluations: 100 -2lnl: 1227.93774
Number of evaluations: 200 -2lnl: 1227.288078
Number of evaluations: 300 -2lnl: 1227.014115
Number of evaluations: 400 -2lnl: 1226.986396
Number of evaluations: 500 -2lnl: 1226.97569
Number of evaluations: 600 -2lnl: 1227.007945
Number of evaluations: 700 -2lnl: 1226.992285
Number of evaluations: 800 -2lnl: 1227.142454
Number of evaluations: 900 -2lnl: 1226.984833
Number of evaluations: 1000 -2lnl: 1227.016204
Number of evaluations: 1100 -2lnl: 1226.985726
Number of evaluations: 1200 -2lnl: 1227.030263
Number of evaluations: 1300 -2lnl: 1226.99256
Number of evaluations: 1400 -2lnl: 1226.97546
Number of evaluations: 100 -2lnl: 1227.911352
Number of evaluations: 200 -2lnl: 1227.036511
Number of evaluations: 300 -2lnl: 1227.646509
Number of evaluations: 400 -2lnl: 1227.025231
Number of evaluations: 500 -2lnl: 1227.074533
Number of evaluations: 600 -2lnl: 1226.990534
Number of evaluations: 700 -2lnl: 1228.152346Phi.habitat.time.p.habitat
##
Number of evaluations: 100 -2lnl: 1227.79551
Number of evaluations: 200 -2lnl: 1227.120457
Number of evaluations: 300 -2lnl: 1226.906674
Number of evaluations: 400 -2lnl: 1226.799329
Number of evaluations: 500 -2lnl: 1226.769009
Number of evaluations: 600 -2lnl: 1226.766748
Number of evaluations: 700 -2lnl: 1226.767176
Number of evaluations: 800 -2lnl: 1226.839261
Number of evaluations: 900 -2lnl: 1226.769448
Number of evaluations: 1000 -2lnl: 1226.849578
Number of evaluations: 1100 -2lnl: 1226.768915
Number of evaluations: 1200 -2lnl: 1226.971117
Number of evaluations: 1300 -2lnl: 1226.793747
Number of evaluations: 1400 -2lnl: 1226.80454
Number of evaluations: 1500 -2lnl: 1226.766644
Number of evaluations: 100 -2lnl: 1227.604133
Number of evaluations: 200 -2lnl: 1226.815512
Number of evaluations: 300 -2lnl: 1227.74426
Number of evaluations: 400 -2lnl: 1226.81178
Number of evaluations: 500 -2lnl: 1226.980183
Number of evaluations: 600 -2lnl: 1226.777627
Number of evaluations: 700 -2lnl: 1228.149466
Number of evaluations: 800 -2lnl: 1226.800166Phi.habitat.time.p.sex
##
Number of evaluations: 100 -2lnl: 1226.29305
Number of evaluations: 200 -2lnl: 1225.546506
Number of evaluations: 300 -2lnl: 1225.331577
Number of evaluations: 400 -2lnl: 1225.233333
Number of evaluations: 500 -2lnl: 1225.22303
Number of evaluations: 600 -2lnl: 1225.221478
Number of evaluations: 700 -2lnl: 1225.222961
Number of evaluations: 800 -2lnl: 1225.255474
Number of evaluations: 900 -2lnl: 1225.242417
Number of evaluations: 1000 -2lnl: 1225.280464
Number of evaluations: 1100 -2lnl: 1225.239853
Number of evaluations: 1200 -2lnl: 1225.252553
Number of evaluations: 1300 -2lnl: 1225.232904
Number of evaluations: 1400 -2lnl: 1225.239373
Number of evaluations: 1500 -2lnl: 1225.22146
Number of evaluations: 100 -2lnl: 1225.899994
Number of evaluations: 200 -2lnl: 1225.267646
Number of evaluations: 300 -2lnl: 1226.220945
Number of evaluations: 400 -2lnl: 1225.270379
Number of evaluations: 500 -2lnl: 1225.447483
Number of evaluations: 600 -2lnl: 1225.232715
Number of evaluations: 700 -2lnl: 1226.486483
Number of evaluations: 800 -2lnl: 1225.258312Phi.habitat.time.p.time
##
Number of evaluations: 100 -2lnl: 1236.170845
Number of evaluations: 200 -2lnl: 1225.998227
Number of evaluations: 300 -2lnl: 1224.510412
Number of evaluations: 400 -2lnl: 1223.789346
Number of evaluations: 500 -2lnl: 1223.614985
Number of evaluations: 600 -2lnl: 1223.525763
Number of evaluations: 700 -2lnl: 1223.408596
Number of evaluations: 800 -2lnl: 1223.362659
Number of evaluations: 900 -2lnl: 1223.314838
Number of evaluations: 1000 -2lnl: 1223.310042
Number of evaluations: 1100 -2lnl: 1223.333229
Number of evaluations: 1200 -2lnl: 1224.496618
Number of evaluations: 1300 -2lnl: 1223.314412
Number of evaluations: 1400 -2lnl: 1223.694594
Number of evaluations: 1500 -2lnl: 1223.327273
Number of evaluations: 1600 -2lnl: 1223.897422
Number of evaluations: 1700 -2lnl: 1223.326004
Number of evaluations: 1800 -2lnl: 1223.548027
Number of evaluations: 1900 -2lnl: 1223.314272
Number of evaluations: 2000 -2lnl: 1223.381878
Number of evaluations: 2100 -2lnl: 1223.32439
Number of evaluations: 2200 -2lnl: 1223.714295
Number of evaluations: 2300 -2lnl: 1223.314922
Number of evaluations: 2400 -2lnl: 1223.331151
Number of evaluations: 2500 -2lnl: 1223.310013
Number of evaluations: 100 -2lnl: 1223.738771
Number of evaluations: 200 -2lnl: 1223.39916
Number of evaluations: 300 -2lnl: 1223.366651
Number of evaluations: 400 -2lnl: 1223.43491
Number of evaluations: 500 -2lnl: 1225.126663
Number of evaluations: 600 -2lnl: 1223.602568
Number of evaluations: 700 -2lnl: 1224.073011
Number of evaluations: 800 -2lnl: 1223.329888
Number of evaluations: 900 -2lnl: 1223.726268
Number of evaluations: 1000 -2lnl: 1223.552908
Number of evaluations: 1100 -2lnl: 1223.425592
Number of evaluations: 1200 -2lnl: 1223.32932
Number of evaluations: 1300 -2lnl: 1224.28456Phi.sex.p.dot
##
Number of evaluations: 100 -2lnl: 1284.628321Phi.sex.p.habitat
##
Number of evaluations: 100 -2lnl: 1284.362121
Number of evaluations: 200 -2lnl: 1284.362121Phi.sex.p.sex
##
Number of evaluations: 100 -2lnl: 1283.680512
Number of evaluations: 200 -2lnl: 1283.680512Phi.sex.p.time
##
Number of evaluations: 100 -2lnl: 1285.058065
Number of evaluations: 200 -2lnl: 1278.590104
Number of evaluations: 300 -2lnl: 1278.287167
Number of evaluations: 400 -2lnl: 1278.550077
Number of evaluations: 500 -2lnl: 1278.274914
Number of evaluations: 600 -2lnl: 1278.278464
Number of evaluations: 100 -2lnl: 1278.487128
Number of evaluations: 200 -2lnl: 1278.281974Phi.sex.time.p.dot
##
Number of evaluations: 100 -2lnl: 1226.568152
Number of evaluations: 200 -2lnl: 1225.990052
Number of evaluations: 300 -2lnl: 1225.526691
Number of evaluations: 400 -2lnl: 1225.520238
Number of evaluations: 500 -2lnl: 1225.441401
Number of evaluations: 600 -2lnl: 1225.457658
Number of evaluations: 700 -2lnl: 1227.229984
Number of evaluations: 800 -2lnl: 1225.823925
Number of evaluations: 900 -2lnl: 1226.666332
Number of evaluations: 1000 -2lnl: 1225.557911
Number of evaluations: 1100 -2lnl: 1226.41834
Number of evaluations: 1200 -2lnl: 1225.546468
Number of evaluations: 1300 -2lnl: 1225.441384
Number of evaluations: 100 -2lnl: 1226.384775
Number of evaluations: 200 -2lnl: 1225.590272
Number of evaluations: 300 -2lnl: 1227.383442
Number of evaluations: 400 -2lnl: 1225.783639
Number of evaluations: 500 -2lnl: 1226.21125
Number of evaluations: 600 -2lnl: 1225.580664
Number of evaluations: 700 -2lnl: 1226.781155Phi.sex.time.p.habitat
##
Number of evaluations: 100 -2lnl: 1226.27436
Number of evaluations: 200 -2lnl: 1225.705522
Number of evaluations: 300 -2lnl: 1225.343191
Number of evaluations: 400 -2lnl: 1225.095631
Number of evaluations: 500 -2lnl: 1225.060522
Number of evaluations: 600 -2lnl: 1225.052203
Number of evaluations: 700 -2lnl: 1225.070844
Number of evaluations: 800 -2lnl: 1226.090349
Number of evaluations: 900 -2lnl: 1225.070343
Number of evaluations: 1000 -2lnl: 1225.434673
Number of evaluations: 1100 -2lnl: 1225.075408
Number of evaluations: 1200 -2lnl: 1225.277906
Number of evaluations: 1300 -2lnl: 1225.063087
Number of evaluations: 1400 -2lnl: 1225.427794
Number of evaluations: 1500 -2lnl: 1225.059431
Number of evaluations: 1600 -2lnl: 1225.05194
Number of evaluations: 100 -2lnl: 1225.877588
Number of evaluations: 200 -2lnl: 1225.472542
Number of evaluations: 300 -2lnl: 1229.615687
Number of evaluations: 400 -2lnl: 1225.117012
Number of evaluations: 500 -2lnl: 1226.496245
Number of evaluations: 600 -2lnl: 1225.122785
Number of evaluations: 700 -2lnl: 1226.623784
Number of evaluations: 800 -2lnl: 1225.094514Phi.sex.time.p.sex
##
Number of evaluations: 100 -2lnl: 1226.094895
Number of evaluations: 200 -2lnl: 1225.21013
Number of evaluations: 300 -2lnl: 1224.774306
Number of evaluations: 400 -2lnl: 1224.607338
Number of evaluations: 500 -2lnl: 1224.548274
Number of evaluations: 600 -2lnl: 1224.538384
Number of evaluations: 700 -2lnl: 1224.787892
Number of evaluations: 800 -2lnl: 1224.598502
Number of evaluations: 900 -2lnl: 1225.403255
Number of evaluations: 1000 -2lnl: 1224.579718
Number of evaluations: 1100 -2lnl: 1224.871416
Number of evaluations: 1200 -2lnl: 1224.548901
Number of evaluations: 1300 -2lnl: 1225.183338
Number of evaluations: 1400 -2lnl: 1224.563466
Number of evaluations: 1500 -2lnl: 1224.538308
Number of evaluations: 1600 -2lnl: 1224.538307
Number of evaluations: 100 -2lnl: 1225.284145
Number of evaluations: 200 -2lnl: 1225.049219
Number of evaluations: 300 -2lnl: 1229.633306
Number of evaluations: 400 -2lnl: 1224.615562
Number of evaluations: 500 -2lnl: 1226.253282
Number of evaluations: 600 -2lnl: 1224.62377
Number of evaluations: 700 -2lnl: 1226.119317
Number of evaluations: 800 -2lnl: 1224.588959Phi.sex.time.p.time
##
Number of evaluations: 100 -2lnl: 1234.690315
Number of evaluations: 200 -2lnl: 1225.62689
Number of evaluations: 300 -2lnl: 1224.314752
Number of evaluations: 400 -2lnl: 1223.075053
Number of evaluations: 500 -2lnl: 1222.724628
Number of evaluations: 600 -2lnl: 1222.508521
Number of evaluations: 700 -2lnl: 1222.332531
Number of evaluations: 800 -2lnl: 1222.210487
Number of evaluations: 900 -2lnl: 1222.133208
Number of evaluations: 1000 -2lnl: 1222.130014
Number of evaluations: 1100 -2lnl: 1222.199949
Number of evaluations: 1200 -2lnl: 1246.351238
Number of evaluations: 1300 -2lnl: 1223.169168
Number of evaluations: 1400 -2lnl: 1239.600774
Number of evaluations: 1500 -2lnl: 1222.436111
Number of evaluations: 1600 -2lnl: 1224.264099
Number of evaluations: 1700 -2lnl: 1222.224401
Number of evaluations: 1800 -2lnl: 1235.88525
Number of evaluations: 1900 -2lnl: 1222.162762
Number of evaluations: 2000 -2lnl: 1222.701721
Number of evaluations: 2100 -2lnl: 1222.194462
Number of evaluations: 2200 -2lnl: 1223.310846
Number of evaluations: 2300 -2lnl: 1222.26185
Number of evaluations: 2400 -2lnl: 1222.129981
Number of evaluations: 2500 -2lnl: 1222.129981
Number of evaluations: 100 -2lnl: 1222.614653
Number of evaluations: 200 -2lnl: 1223.472304
Number of evaluations: 300 -2lnl: 1225.226947
Number of evaluations: 400 -2lnl: 1222.614453
Number of evaluations: 500 -2lnl: 1224.253244
Number of evaluations: 600 -2lnl: 1222.453988
Number of evaluations: 700 -2lnl: 1226.374387
Number of evaluations: 800 -2lnl: 1222.358573
Number of evaluations: 900 -2lnl: 1225.989839
Number of evaluations: 1000 -2lnl: 1222.382887
Number of evaluations: 1100 -2lnl: 1222.225911
Number of evaluations: 1200 -2lnl: 1222.161356
Number of evaluations: 1300 -2lnl: 1223.262178Phi.time.p.dot
##
Number of evaluations: 100 -2lnl: 1242.648575
Number of evaluations: 200 -2lnl: 1242.487127
Number of evaluations: 300 -2lnl: 1245.445377
Number of evaluations: 400 -2lnl: 1242.418534
Number of evaluations: 100 -2lnl: 1243.038192
Number of evaluations: 200 -2lnl: 1242.446706Phi.time.p.habitat
##
Number of evaluations: 100 -2lnl: 1242.360719
Number of evaluations: 200 -2lnl: 1241.978159
Number of evaluations: 300 -2lnl: 1243.089518
Number of evaluations: 400 -2lnl: 1241.989515
Number of evaluations: 500 -2lnl: 1242.041195
Number of evaluations: 100 -2lnl: 1243.15084
Number of evaluations: 200 -2lnl: 1242.099554Phi.time.p.sex
##
Number of evaluations: 100 -2lnl: 1240.930825
Number of evaluations: 200 -2lnl: 1240.477856
Number of evaluations: 300 -2lnl: 1242.083998
Number of evaluations: 400 -2lnl: 1240.508214
Number of evaluations: 500 -2lnl: 1241.395697
Number of evaluations: 100 -2lnl: 1241.640046
Number of evaluations: 200 -2lnl: 1240.590849Phi.time.p.time
##
Number of evaluations: 100 -2lnl: 1244.670707
Number of evaluations: 200 -2lnl: 1239.725494
Number of evaluations: 300 -2lnl: 1239.319519
Number of evaluations: 400 -2lnl: 1239.303701
Number of evaluations: 500 -2lnl: 1239.350357
Number of evaluations: 600 -2lnl: 1239.417367
Number of evaluations: 700 -2lnl: 1243.575689
Number of evaluations: 800 -2lnl: 1239.311438
Number of evaluations: 900 -2lnl: 1241.176613
Number of evaluations: 1000 -2lnl: 1239.329909
Number of evaluations: 1100 -2lnl: 1239.303675
Number of evaluations: 100 -2lnl: 1247.504101
Number of evaluations: 200 -2lnl: 1239.375924
Number of evaluations: 300 -2lnl: 1240.074042
Number of evaluations: 400 -2lnl: 1239.422373
Number of evaluations: 500 -2lnl: 1239.433946
Number of evaluations: 600 -2lnl: 1239.305696
Look at model results
summary(starlings_cjs_models)## Length Class Mode
## Phi.dot.p.dot 5 crm list
## Phi.dot.p.habitat 5 crm list
## Phi.dot.p.sex 5 crm list
## Phi.dot.p.time 5 crm list
## Phi.habitat.time.p.dot 5 crm list
## Phi.habitat.time.p.habitat 5 crm list
## Phi.habitat.time.p.sex 5 crm list
## Phi.habitat.time.p.time 5 crm list
## Phi.sex.p.dot 5 crm list
## Phi.sex.p.habitat 5 crm list
## Phi.sex.p.sex 5 crm list
## Phi.sex.p.time 5 crm list
## Phi.sex.time.p.dot 5 crm list
## Phi.sex.time.p.habitat 5 crm list
## Phi.sex.time.p.sex 5 crm list
## Phi.sex.time.p.time 5 crm list
## Phi.time.p.dot 5 crm list
## Phi.time.p.habitat 5 crm list
## Phi.time.p.sex 5 crm list
## Phi.time.p.time 5 crm list
## model.table 7 data.frame list
Time-dependent variables
Heat_wave <- matrix(rep(c(0, 1, 1, 0, 0, 0), each = nrow(starlings)), ncol =
6)
colnames(Heat_wave) = paste("Heat_wave", 1:6, sep = "")
starlings_heatwave <- cbind(starlings, Heat_wave)
#generate some random weight
starlings_heatwave$weight <- round(rnorm(nrow(starlings), 80, 3))
#prepare data
starlings.proc <- process.data(starlings_heatwave)
design.Phi <-
list(static = c("weight", "sex"),
time.varying = c("Heat_wave"))
design.p <- list(static = c("sex"))
design.parameters <- list(Phi = design.Phi, p = design.p)
m3_design <- make.design.data(data = starlings.proc,
parameters = design.parameters)
names(m3_design$Phi)## [1] "id" "occ" "time" "cohort" "age" "Heat_wave"
## [7] "weight" "sex" "Time" "Cohort" "Age" "order"
names(m3_design$p)## [1] "id" "occ" "time" "cohort" "age" "sex" "Time" "Cohort"
## [9] "Age" "fix" "order"
Phi.heat.sex <- list(formula = ~ Heat_wave * sex)
p.sex <- list(formula = ~ sex)
m3 <- crm(
starlings.proc,
m3_design,
hessian = TRUE,
model.parameters = list(Phi = Phi.heat.sex,
p = p.sex)
)##
Number of evaluations: 100 -2lnl: 1244.121204
Number of evaluations: 200 -2lnl: 1243.982328
Number of evaluations: 300 -2lnl: 1244.263347
Number of evaluations: 100 -2lnl: 1245.196681
predict(m3)## $Phi
## Heat_wave sex occ estimate se lcl ucl
## 1 0 Female 6 0.5888999 0.03365102 0.5217320 0.6529102
## 2 0 Male 6 0.6029780 0.03050136 0.5419426 0.6609683
## 3 1 Female 3 0.2874165 0.03769815 0.2194742 0.3665152
## 4 1 Male 3 0.4682230 0.04526839 0.3813970 0.5570169
##
## $p
## sex occ estimate se lcl ucl
## 1 Female 7 0.8722827 0.03817459 0.7772436 0.9304042
## 2 Male 7 0.9191637 0.02691341 0.8482694 0.9585520
More complex models & predicting
Phi.weight.heat <-list(formula=~weight*Heat_wave)
p.sex <-list(formula=~sex)
m4 <- crm(
starlings.proc,
m3_design,
hessian = TRUE,
model.parameters = list(Phi = Phi.weight.heat,
p = p.sex)
)##
Number of evaluations: 100 -2lnl: 1252.677436
Number of evaluations: 200 -2lnl: 1255.89512
Number of evaluations: 300 -2lnl: 1252.968609
Number of evaluations: 400 -2lnl: 1252.676557
Number of evaluations: 100 -2lnl: 1257.437678
#predict model output
predict(m4)## $Phi
## weight Heat_wave occ estimate se lcl ucl
## 1 79 0 6 0.6003154 0.02366211 0.5531702 0.6456726
## 2 78 0 6 0.6030211 0.02649719 0.5501114 0.6536262
## 3 76 0 6 0.6084135 0.03577190 0.5365250 0.6758874
## 4 87 0 6 0.5784621 0.05581375 0.4670070 0.6824580
## 5 84 0 6 0.5866981 0.03700753 0.5127970 0.6568898
## 6 77 0 6 0.6057205 0.03070027 0.5442351 0.6640317
## 7 82 0 6 0.5921622 0.02713700 0.5380946 0.6440871
## 8 80 0 6 0.5976036 0.02274455 0.5523348 0.6412679
## 9 81 0 6 0.5948858 0.02399419 0.5471243 0.6409161
## 10 85 0 6 0.5839579 0.04294503 0.4981462 0.6649651
## 11 75 0 6 0.6110998 0.04137545 0.5276484 0.6885113
## 12 83 0 6 0.5894329 0.03163288 0.5263314 0.6497225
## 13 86 0 6 0.5812125 0.04925269 0.4827940 0.6735658
## 14 73 0 6 0.6164522 0.05343803 0.5078826 0.7145325
## 15 74 0 4 0.6137795 0.04730548 0.5180168 0.7014801
## 16 90 0 4 0.5701823 0.07641024 0.4186053 0.7096513
## 17 80 1 3 0.3691765 0.02984352 0.3128306 0.4293305
## 18 86 1 3 0.4125508 0.06882940 0.2869833 0.5506308
## 19 78 1 3 0.3551391 0.03739509 0.2856581 0.4313161
## 20 77 1 3 0.3482111 0.04407812 0.2674567 0.4387446
## 21 81 1 3 0.3762808 0.03081215 0.3180633 0.4383054
## 22 79 1 3 0.3621284 0.03230885 0.3014750 0.4275162
## 23 84 1 3 0.3979045 0.04992160 0.3052082 0.4985524
## 24 83 1 3 0.3906476 0.04173233 0.3125535 0.4747780
## 25 82 1 3 0.3834388 0.03507171 0.3174010 0.4540764
## 26 76 1 3 0.3413467 0.05162750 0.2483609 0.4483786
## 27 87 1 3 0.4199341 0.07901220 0.2771182 0.5775476
## 28 85 1 3 0.4052065 0.05906622 0.2964706 0.5241130
## 29 75 1 2 0.3345482 0.05961159 0.2292578 0.4593738
## 30 88 0 6 0.5757068 0.06255564 0.4509656 0.6914950
## 31 72 0 5 0.6191180 0.05969651 0.4974060 0.7275026
## 32 74 1 3 0.3278177 0.06778289 0.2106540 0.4712427
## 33 72 1 2 0.3145687 0.08415756 0.1759748 0.4965423
##
## $p
## sex occ estimate se lcl ucl
## 1 Female 7 0.8568973 0.03983628 0.7600623 0.9188251
## 2 Male 7 0.9240681 0.02511068 0.8578462 0.9608486
new_starling <- expand.grid(sex=c("Female","Male"),weight=60:80,Heat_wave1=0,
Heat_wave2=1,Heat_wave3=1,Heat_wave4=0,Heat_wave5=0,Heat_wave6=0)
prediction <- predict(m4,newdata=new_starling,se=TRUE)
prediction$Phi$Heat_wave = factor(prediction$Phi$Heat_wave, labels = c("No_Heatwave", "Heatwave"))
ggplot(prediction$Phi, aes(weight, estimate, ymin = lcl, ymax = ucl)) +
geom_errorbar(width = 0.2) +
geom_point() +
geom_line() +
xlab("\nWeight") + ylab("Survival\n") + facet_grid(Heat_wave ~ .)Exercise
Some exercises:
Exercise 6.1
The data set starlings contains one column (“infected”) that represents an infection with an incurable disease (i.e. a static variable comparable to sex). Formulate hypothesis associated with this infection status (write them down!) & try to build simple CMR models. Think about why infection could matter for capture probability or survival. There are some NAs in the data set, maybe try ?na-omit() and remove some of the observations without data from all analyses.Discuss the results in one sentence
Exercise 6.2
In our previous model (m3), we saw that sex & heatwave played an important role in individual survival. Formulate a model that additionally accounts for the habitat type (hint: * or :)
Exercise 6.3 (advanced)
This exercise will require some data-wrangling, the data set is not yet prepared: In the file locations, you will find measured temperatures for each year.
(Hint: You might need dplyr::left_join() and dplyr::rename() to get all data in the right format and bind the rows according to the habitat type.)
Test how temperature affects survival, and how sex affects detection probability. Then predict the survival probabilities for extreme weather scenarios (temperatures of 30,33,35,38,41,44 degree Celsius)
Solution 6.1
There are many solutions, one is:
Hypothesis: infection affects survival & prediction
Prediction: Lower survival AND detection probability in infected individuals. Here, we hypothesize that infected animals move less, therefore are less likely to be captured. Additionally, infected animals suffer from disease, thus survival is decreased
package.list=c("here",
"marked",
"tidyverse")
for (package in package.list) {
if (!require(package, character.only=T, quietly=T)) {
install.packages(package)
library(package, character.only=T)
}
}
starlings <-
read_csv(here("data", "data_berlin", "animal_data", "starlings.csv")) %>%
mutate(
sex = as.factor(sex),
habitat = as.factor(habitat),
infected = as.factor(infected)
) %>%
rename(ID = "...1")%>%
as.data.frame()
#load environmental data on capturing
locations <-
read_csv(here(
"data",
"data_berlin",
"animal_data",
"starlings_capture_location_data.csv"
))
temperatures <- locations %>%
dplyr::select(-c("Lon", "Lat")) %>%
rename(habitat = Habitat)
m_infection<- process.data(data = na.omit(starlings),groups = c("infected"))
m_infetion_design <- make.design.data(data = m_infection)
phi_infection <- list(formula = ~ infected)
p_infection <- list(formula = ~ infected)
cjs_infection <- crm(
m_infection,
m_infetion_design,
model.parameters = list(Phi = phi_infection,
p = p_infection),
accumulate = FALSE
)##
Number of evaluations: 100 -2lnl: 861.9151509
cjs_infection##
## crm Model Summary
##
## Npar : 4
## -2lnL: 861.9147
## AIC : 869.9147
##
## Beta
## Estimate
## Phi.(Intercept) -0.03261002
## Phi.infectedYes 0.20133675
## p.(Intercept) 2.68053887
## p.infectedYes -0.76975772
predict(object = cjs_infection)## $Phi
## infected occ estimate
## 1 Yes 6 0.5420819
## 2 No 6 0.4918482
##
## $p
## infected occ estimate
## 1 Yes 7 0.8711069
## 2 No 7 0.9358685
Solution 6.2
#prepare data
Heat_wave <- matrix(rep(c(0, 1, 1, 0, 0, 0), each = nrow(starlings)), ncol =
6)
colnames(Heat_wave) = paste("Heat_wave", 1:6, sep = "")
starlings_heatwave <- cbind(starlings, Heat_wave)
#generate some random weight
starlings_heatwave$weight <- round(rnorm(nrow(starlings), 80, 3))
starlings.proc <- process.data(starlings_heatwave)
design.Phi <-
list(static = c("weight", "sex","habitat"),
time.varying = c("Heat_wave"))
design.p <- list(static = c("sex"))
design.parameters <- list(Phi = design.Phi, p = design.p)
m6_design <- make.design.data(data = starlings.proc,
parameters = design.parameters)
names(m6_design$Phi)## [1] "id" "occ" "time" "cohort" "age" "Heat_wave"
## [7] "weight" "sex" "habitat" "Time" "Cohort" "Age"
## [13] "order"
names(m6_design$p)## [1] "id" "occ" "time" "cohort" "age" "sex" "Time" "Cohort"
## [9] "Age" "fix" "order"
Phi.heat.sex.habitat <- list(formula = ~ Heat_wave * sex * habitat)
p.sex <- list(formula = ~ sex)
m6 <- crm(
starlings.proc,
m6_design,
hessian = TRUE,
model.parameters = list(Phi = Phi.heat.sex.habitat,
p = p.sex)
)##
Number of evaluations: 100 -2lnl: 1224.55332
Number of evaluations: 200 -2lnl: 1224.114866
Number of evaluations: 300 -2lnl: 1224.07174
Number of evaluations: 400 -2lnl: 1224.076663
Number of evaluations: 500 -2lnl: 1224.093993
Number of evaluations: 600 -2lnl: 1224.305982
Number of evaluations: 700 -2lnl: 1224.072706
Number of evaluations: 800 -2lnl: 1224.071643
Number of evaluations: 100 -2lnl: 1224.265117
Number of evaluations: 200 -2lnl: 1224.137044
Number of evaluations: 300 -2lnl: 1224.878393
Number of evaluations: 400 -2lnl: 1224.083587
predict(m6)## $Phi
## Heat_wave sex habitat occ estimate se lcl ucl
## 1 0 Female rural 6 0.5359784 0.04604001 0.44554937 0.6241018
## 2 0 Male rural 6 0.5956507 0.04259539 0.51018411 0.6756843
## 3 1 Female rural 3 0.4429930 0.05930527 0.33178993 0.5602196
## 4 1 Male rural 3 0.4807007 0.06379440 0.35936302 0.6043584
## 5 0 Female urban 6 0.6391793 0.04551414 0.54611911 0.7228425
## 6 0 Male urban 6 0.6103285 0.04248132 0.52462849 0.6897163
## 7 1 Female urban 3 0.1365968 0.04010417 0.07513682 0.2355269
## 8 1 Male urban 3 0.4555735 0.06396440 0.33544774 0.5811012
##
## $p
## sex occ estimate se lcl ucl
## 1 Female 7 0.8739003 0.03775411 0.7797819 0.9313357
## 2 Male 7 0.9191464 0.02691835 0.8482405 0.9585424
Solution 6.3
starlings_temperatur <- starlings %>%
left_join(temperatures, by = "habitat")%>%
rename(temperature1=temperature_2016,
temperature2=temperature_2017,
temperature3=temperature_2018,
temperature4=temperature_2019,
temperature5=temperature_2020,
temperature6=temperature_2021,
temperature7=temperature_2022)
#prepare data
starlings.proc.temp <- process.data(starlings_temperatur)
design.Phi <-
list(static = c("sex"),
time.varying = c("temperature"))
design.p <- list(static = c("sex"))
design.parameters <- list(Phi = design.Phi, p = design.p)
m7_design <- make.design.data(data = starlings.proc.temp,
parameters = design.parameters)
names(m7_design$Phi)## [1] "id" "occ" "time" "cohort" "age"
## [6] "temperature" "sex" "Time" "Cohort" "Age"
## [11] "order"
names(m7_design$p)## [1] "id" "occ" "time" "cohort" "age" "sex" "Time" "Cohort"
## [9] "Age" "fix" "order"
Phi.heat.sex <- list(formula = ~ temperature)
p.sex <- list(formula = ~ sex)
m7 <- crm(
starlings.proc.temp,
m7_design,
hessian = TRUE,
model.parameters = list(Phi = Phi.heat.sex,
p = p.sex)
)##
Number of evaluations: 100 -2lnl: 1249.569476
Number of evaluations: 200 -2lnl: 1242.783021
predict(m7)## $Phi
## temperature occ estimate se lcl ucl
## 1 29.0 6 0.5900881 0.02104224 0.5483047 0.6306103
## 2 27.0 4 0.6457556 0.02513997 0.5950862 0.6933517
## 3 31.0 3 0.5320125 0.01892171 0.4948173 0.5688552
## 4 35.0 2 0.4148422 0.02469059 0.3674108 0.4639063
## 5 29.5 1 0.5757385 0.02023851 0.5356662 0.6148391
## 6 28.0 5 0.6183101 0.02298216 0.5723670 0.6622301
## 7 30.0 4 0.5612600 0.01959489 0.5225627 0.5992263
## 8 36.0 3 0.3865030 0.02722904 0.3346766 0.4410349
## 9 43.0 2 0.2161236 0.03866961 0.1498497 0.3013197
##
## $p
## sex occ estimate se lcl ucl
## 1 Female 7 0.8699031 0.03750525 0.7773961 0.9275499
## 2 Male 7 0.9266009 0.02440500 0.8620348 0.9622732
new_starling2 <- expand.grid(sex=c("Female","Male"),temperature1=30,
temperature2=33,temperature3=35,temperature4=38,temperature5=41,temperature6=44,weight=60:80)
prediction <- predict(m7,newdata=new_starling2,se=TRUE)
prediction$Phi$temperature = factor(prediction$Phi$temperature)
ggplot(prediction$Phi, aes(temperature, estimate, ymin = lcl, ymax = ucl)) +
geom_errorbar(width = 0.2) +
geom_point() +
geom_line() +
xlab("\nTemperature") + ylab("Survival probability\n")Session Info
Sys.time()## [1] "2023-03-16 19:09:26 CET"
#git2r::repository() ## uncomment if you are using GitHub
sessionInfo()## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8
## [2] LC_CTYPE=English_United Kingdom.utf8
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=C
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] tmap_3.3-3 sf_1.0-9 forcats_0.5.2 stringr_1.5.0
## [5] dplyr_1.0.10 purrr_1.0.1 readr_2.1.3 tidyr_1.3.0
## [9] tibble_3.1.8 ggplot2_3.4.0 tidyverse_1.3.2 skimr_2.1.5
## [13] marked_1.2.6 lme4_1.1-31 Matrix_1.5-3 here_1.0.1
##
## loaded via a namespace (and not attached):
## [1] leafem_0.2.0 googledrive_2.0.0 minqa_1.2.5
## [4] colorspace_2.1-0 ellipsis_0.3.2 class_7.3-21
## [7] leaflet_2.1.1 rprojroot_2.0.3 base64enc_0.1-3
## [10] fs_1.6.0 dichromat_2.0-0.1 rstudioapi_0.14
## [13] proxy_0.4-27 farver_2.1.1 bit64_4.0.5
## [16] fansi_1.0.4 lubridate_1.9.1 xml2_1.3.3
## [19] codetools_0.2-18 splines_4.2.2 cachem_1.0.6
## [22] knitr_1.41 jsonlite_1.8.4 nloptr_2.0.3
## [25] tmaptools_3.1-1 broom_1.0.2 dbplyr_2.3.0
## [28] png_0.1-8 compiler_4.2.2 httr_1.4.4
## [31] backports_1.4.1 assertthat_0.2.1 fastmap_1.1.0
## [34] gargle_1.2.1 cli_3.6.0 s2_1.1.2
## [37] leaflet.providers_1.9.0 htmltools_0.5.4 tools_4.2.2
## [40] coda_0.19-4 gtable_0.3.1 glue_1.6.2
## [43] wk_0.7.1 Rcpp_1.0.10 raster_3.6-14
## [46] cellranger_1.1.0 jquerylib_0.1.4 vctrs_0.5.2
## [49] nlme_3.1-161 leafsync_0.1.0 crosstalk_1.2.0
## [52] lwgeom_0.2-11 xfun_0.36 rvest_1.0.3
## [55] timechange_0.2.0 lifecycle_1.0.3 XML_3.99-0.13
## [58] googlesheets4_1.0.1 terra_1.7-3 MASS_7.3-58.2
## [61] scales_1.2.1 vroom_1.6.1 ragg_1.2.5
## [64] hms_1.1.2 expm_0.999-7 TMB_1.9.2
## [67] RColorBrewer_1.1-3 yaml_2.3.7 sass_0.4.5
## [70] stringi_1.7.12 highr_0.10 e1071_1.7-12
## [73] optimx_2022-4.30 boot_1.3-28.1 truncnorm_1.0-8
## [76] R2admb_0.7.16.3 repr_1.1.5 rlang_1.0.6
## [79] pkgconfig_2.0.3 systemfonts_1.0.4 evaluate_0.20
## [82] lattice_0.20-45 labeling_0.4.2 htmlwidgets_1.6.1
## [85] bit_4.0.5 tidyselect_1.2.0 magrittr_2.0.3
## [88] bookdown_0.32 R6_2.5.1 generics_0.1.3
## [91] DBI_1.1.3 pillar_1.8.1 haven_2.5.1
## [94] withr_2.5.0 units_0.8-1 stars_0.6-0
## [97] sp_1.6-0 abind_1.4-5 modelr_0.1.10
## [100] crayon_1.5.2 KernSmooth_2.23-20 utf8_1.2.2
## [103] tzdb_0.3.0 rmarkdown_2.20 grid_4.2.2
## [106] readxl_1.4.1 data.table_1.14.6 rmdformats_1.0.4
## [109] reprex_2.0.2 digest_0.6.31 classInt_0.4-8
## [112] numDeriv_2016.8-1.1 textshaping_0.3.6 munsell_0.5.0
## [115] viridisLite_0.4.1 bslib_0.4.2